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ABSTRACT 

Pluto's recently discovered minor moons, Nix and Hydra, have almost circular orbits, and are nearly 
coplanar with Charon, Pluto's major moon. This is surprising because tidal interactions with Pluto 
arc too weak to damp their eccentricities. Wc consider an alternative possibility: that Nix and Hydra 
circularize their orbits by exciting Charon's eccentricity via secular interactions, and Charon in turn 
damps its own eccentricity by tidal interaction with Pluto. The timescale for this process can be less 
than the age of the Solar System, for plausible tidal parameters and moon masses. However, as we 
show numerically and analytically, the effects of the 2:1 and 3:1 resonant forcing terms between Nix 
and Charon complicate this picture. In the presence of Charon's tidal damping, the 2:1 term forces 
Nix to migrate outward and the 3:1 term changes the eccentricity damping rate, sometimes leading 
to eccentricity growth. We conclude that this mechanism probably does not explain Nix and Hydra's 
current orbits. Instead, we suggest that they were formed in-situ with low eccentricities. 

We also show that an upper limit on Nix's migration speed sets a lower limit on Pluto-Charon's 
tidal circularization timescale of > 10^ yrs. Moreover, Hydra's observed proper eccentricity may be 
explained by the 3:2 forcing by Nix. 
Subject headings: 



1. INTRODUCTION 

Weaver et al. (2006) discovered two small moons or- 
biting Pluto. These moons — Nix and Hydra — are much 

less massive than Pluto's major moon Charon, whose 
mass relative to Pluto's is Mc/Mp ~ 0.12 (Buie et al. 
2006), whereas Nix and Hydra's masses are M^/Mp = 
(4±4) X 10-5 and Mh/Mp = (2±4) x 10-^ based on a 4- 
body fit to the observed positions (Tholen et al. 2007b). ^ 
Weaver et al. (2006) also found that (a) the orbits of Nix 
and Hydra are nearly circular and coplanar with Charon; 
and (b) the period ratios of Charon:Nix:Hydra are nearly 
1:4:6. By analyzing earlier data taken over a year-long 
interval and fitting the data to Keplerian orbits around 
a point mass, Buie et al. (2006) showed that (a) Nix 
and Charon's eccentricities are consistent with zero, with 
CN = 0.0023(21) and ec = 0.00000(7) where brackets de- 
note 1-cr errors in the trailing digits, but Hydra's formally 
is not, with e^r = 0.0052(11); and (b) the period ratio 
of Charon:Nbc:Hydra is 1:3.8915(2):5.9817(2). Recently 
Tholen et al. (2007b) performed a full four-body fit to the 
moons' positions, and thereby obtained slightly different 
orbital parameters. We discuss some of their results in 
§3. 

Hopefully, Nix and Hydra's remarkable orbital prop- 
erties can teach us something about how they formed, 

perhaps shedding light on the formation of Kuiper belt 
objects in general. The near circularity and coplanarity is 
surprising. If Nix and Hydra formed in the collision that 
produced Charon, they likely had high eccentricities and 
inclinations just after formation. Nix and Hydra cannot 
damp their eccentricities through tidal interactions with 
Pluto — the corresponding damping times, for typical pa- 
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^ Although these masses are consistent with zero, the bright- 
ness of the moons implies that Mi^,Mh ~ 1 — 50 x W~^Mp for 
reasonable densities and albedoes (Weaver et al. 2006). 



rameters, are at least 10^ times longer than the age of 
the Solar System (also see Stern et al. 2006).^ 

Ward & Canup (2006) propose a scenario that not only 
accounts for Nix and Hydra's low eccentricities and incli- 
nations, but also accounts for their near-resonant orbits. 
In their "forced resonant migration" scenario. Nix and 
Hydra formed in the impact that produced Charon. Im- 
mediately after impact, all three moons were nnich closer 
to Pluto than they are today, and Nix and Hydra's ec- 
centricities were damped by a disk of post-impact debris. 
Subsequently, Charon migrated outwards by raising tides 
on Pluto, and it forced Nix and Hydra to migrate as well 
because they were trapped in its 4:1 and 6:1 corotation 
resonances. Just before Charon stopped migrating, it left 
Nix and Hydra at their current positions. Although this 
is a neat scenario, we show in a forthcoming paper (Lith- 
wich & Wu, in preparation) that it cannot work. The dif- 
ficulty is that Charon cannot simultaneously force both 
Hydra and Nix to migrate. To transport Nix, Charon's 
eccentricity must satisfy ec < 0.024; otherwise, the 4:1 
corotation resonance is destroyed by resonance overlap, 
as we prove with numerical simulations. But to trans- 
port Hydra, it is required that ec > 0.1; otherwise, the 
width of the 6:1 resonance is so small that the time for 
the resonance to sweep across Hydra is shorter than the 
libration time within the resonance. Because these two 
constraints are incompatible, the Ward & Canup (2006) 
scenario cannot work. 

But Nix and Hydra can damp their eccentricities in a 
different way that has not been considered previously: 
they can transfer their eccentricity to Charon via secular 
interactions, and since Charon damps its own eccentricity 
relatively quickly through its own tidal interactions with 
Pluto, this process can damp eccentricities faster than 

^ This is assuming monolithic bodies. For Nix, the e-folding time 
can be brought to within a few Gyrs if it is a strcngthlcss rubble 
pile and has both tidal Love number and tidal Q factor of order 
unity. 
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the direct interaction of Nix or Hydra with Pluto. We 
shall show below that the timescale for this process for 
Nix is (eq. [24]) 

where tc is Charon's eccentricity damping time. The 
value of Tc is somewhat uncertain. If one assumes that 
Pluto and Charon arc consolidated ice spheres with ma- 
terial strength of 4 x 10^'^erg/cm"^, then tc ~ 5(Q/100) 
Myr (Goldrcich 1963; Dobrovolskis et al. 1997), where 
Q is either Pluto's or Charon's tidal damping parame- 
ter. (Both Pluto and Charon contribute comparably to 
Tc-) Taking rc = 5 Myr, equation (1) implies that the 
damping timescale is less than the age of the Solar Sys- 
tem if M^/Mp > 3 X 10^^. Given the uncertainties in 
Tc — for example, Q could be <C 100, or Charon might be 
a rubble pile instead of a consolidated ice sphere, which 
would increase its tidal Love number and decrease tc it 
is quite possible that the timescale of equation (1) is less 
than the age of the Solar System. 

One is then tempted to conclude that Nix and Hydra 
could have formed with relatively large eccentricities, and 
that these were damped away as described above. How- 
ever, we shall show below that, for Nix, the 2:1 and 3:1 
resonant forcing terms with Charon severely complicate 
this picture. 

2. PLUTO, CHARON AND NIX 
2.1. N-body Simulation 

Figure 1 shows the result of an N-body simulation 

of Pluto, Charon, and Nix, with tidal damping acting 
on Pluto and Charon. The N-body integration is per- 
formed with the SWIFT package (Levison & Duncan 
1994), using the hierarchical Jacobi symplectic integra- 
tor of Beust (2003), and supplemented with an algo- 
rithm for tidal damping. Nix and Charon have masses 
Mn = 1.5 X IQ-^Mp and Mc = O.lMp, respectively. To 
speed up the simulation, we set the circularization time 
of the Pluto-Charon binary to rc = 100 yrs. Although 
this is shorter than the true tc by many orders of mag- 
nitude, we shall show that all relevant timescales scale 
with Tc- Hence for more realistic values of rc, one need 
only rescale the bottom time axis by the factor tc /WO 
yrs. 

The points in the top panel show Charon and Nix's 

instantaneous eccentricities in Jacobi coordinates, and 
the lines show their proper eccentricities. To be more 
precise, the lines are the time-average of one component 
of the eccentricity vectors, averaged over 100 days, which 
is a crude way to remove the short-term noise from the 
eccentricities. We start Charon on a circular orbit at its 
current location (with a period of 6.4 days), and Nix on 
an orbit with an instantaneous eccentricity 0.05 and with 
dN/dc ~ 2.2, where a is the semi-major axis. At early 
times. Nix's proper eccentricity decays. The top time 
axis is in units of the damping time of the slowly damped 
secular mode (eq. [24], which is equivalent to eq. [1]). 
The figure shows that Nix's eccentricity initially damps 
on a timescale much faster than equation (1). But at 
later times. Nix's eccentricity evolves in a complicated 
manner: it increases for a while before decaying again. 

The bottom panel of Figure 1 shows the ratio of semi- 
major axes. Charon's semimajor axis ac (not shown) 
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Fig. 1. — N-body simulation of Pluto, Charon, and Nix with 
tidal damping on Pluto and Charon. The simulation had tidal 
circularization time tc = lOOyrs; for more realistic values of tq, 
the bottom time axis should be rescaled by rc/lOOyrs. The top 
time axis is in units of the damping time of the slowly damped 
secular mode (§2.3). In the top panel, the points depict Charon 
and Nix's instantaneous eccentricities in Jacobi coordinates, and 
the solid lines depict their proper eccentricities. The lower panel 
displays the semi-major axis ratio between Nbc and Charon, with 
the arrow indicating the nominal position of the 4 : 1 resonance. 

remains virtually constant. Nix's semimajor axis slowly 

increases throughout the simulation, and crosses through 
the nominal position of the 4: 1 resonance (marked by the 
arrow). 

2.2. Simplified Model 

In the remainder of this section, we explain the compli- 
cated behavior shown in Figure 1. In this subsection we 
introduce a simplified model, keeping only a small num- 
ber of terms in the disturbing function, and show that 
a numerical integration of the simplified model qualita- 
tively reproduces the behavior of Figure 1. In §§2.3- 
2.5, we isolate the individual terms within the simplified 
model, and explain their effects. The reader uninter- 
ested in the technical details may skip to a summary of 
our findings in §2.6. 

We adopt Jacobi orbital elements, in which Nix's el- 
ements are relative to the center-of-mass of the Pluto- 
Charon binary, and Charon's elements are relative to 
Pluto: 

{a,,A,,ej,Wj} , j G {C,N} (2) 

where subscript C is for Charon, or to be more accurate 
the Pluto-Charon binary, and A'' is for Nix. Standard 
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treatments of the planetary equations adopt heliocen- 
tric (equivalent here to Plutocentric) elements (Murray & 
Dermott 1999). But the present problem is more suited 
to Jacobi elements. We show in the Appendix that us- 
ing Jacobi elements changes the standard treatment in 
two ways. First, instead of the usual indirect term in 
the disturbing function, there is a different correction to 
the direct term, given by equation (A23) to first order in 
the ratio of Charon's to Pluto's mass {Mc/Mp ~ 0.1). 
And second, the appropriate masses must be used; for 
example, it is the reduced mass of the Pluto-Charon bi- 
nary that enters into equation of motion. We ignore this 
second change in the body of the paper because we only 
seek an accuracy of ~ Mc/Mp ~ 10%.^ 

We model the gravitational forces on Nix and Charon 
with the truncated Hamiltonian 



H = Hunp.C + Hunp.N + Hscc,CN + H2:1.CN + H^-^i^CN ■ 

(3) 

These terms and their associated coefficients are defined 
in Tables 1 and 2, respectively, where Mj are the masses, 
fic = Mc/{Mc + Mp), 

Zj ^ e,e^^^ (4) 

are the complex eccentricities, and we choose units so 
that 

C?(Mp + Mc) = 1 . (5) 

The equations of motion are given by Hamilton's equa- 
tions^ 



dt 

dttj 
~dF 

dt 



2i dH 

I 

2^ dH 



Mj^dz* 



dH 

Mj daj 



ie{C,iV} 



(6) 
(7) 
(8) 



For dXj/dt, it suffices in this paper to consider only the 
contribution of the unperturbed energies H^upj, i-e. 

dXn 

—r- = n-- 
dt 
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where 



_ -3/2 



Tii = a 



(9) 
(10) 



is the mean motion. 

Both the 2:1 and the 3:1 resonant forcing terms in the 
Hamiltonian play important roles even when Nix and 
Charon are not particularly close to the nominal loca- 
tions of those resonances. 

^ Adopting Plutocentric elements would give corrections relative 
to the Jacobi elements that are significantly larger than Mc/Mp. 
To appreciate this, consider what would happen if the Pluto- 
Charon binary was very tight. Then Pluto's large reflex velocity 
relative to the binary's center-of-mass would enter into Nix's Plu- 
tocentric elements, even though Nix would be weakly perturbed 
by the non-pointlike nature of the binary. Using Jacobi elements 
avoids this effect. 

® The exact equations arc given by equations (7)-{8) and, in place 
of equation (6), dCj/dt = -idH/dC, where Cj = i^j^y^^i^ - 



— e?)^/^e'^5 is a complex canonical variable (Ogilvie 2007). 

In converting to equation (6), we set C,j = {Mj^faj/2)^/'^Zj, valid 
to leading order in Cj, and assume that aj is constant, because 
corrections caused by varying aj axe higher order in eccentricity. 
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Fig. 2. — Numerical integration of the simplified model. The 
curves show the result from integrating equations (6), (7), and (9) 
with tidal damping (eq. [11]), and with the Hamiltonian given by 
eq. (3). The masses of Pluto, Charon, and Nix and the tidal damp- 
ing rate of Charon are all the same as in the N-body simulation 
of Figure 1. The top panel shows the eccentricities after subtract- 
ing off the complex eccentricities forced by the 2:1 resonance (eqs. 
[31]-[32]). The evolution seen in this figure is qualitatively similar 
to that seen in Figure 1: the eccentricities first decay, then rise, 
then decay again. And Nix is pushed out by Charon. 

We model the tidal damping of Charon's eccentricity 
by adding the term 



dzc 
dt 



tide 



-iczc (11) 



to the equation for zc, where 



ic = I/tc (12) 

is Charon's eccentricity damping rate in the absence of 
Nix. 

Figtire 2 shows the result from numerically integrating 
these equations of motion, where the gi coefficients in 
the Hamiltonian are evaluated at the 4:1 location (Table 
2). The evolution in this simplified model is qualitatively 
similar to that seen in the N-body integration of Figure 
1. The eccentricities initially decay, then rise, and finally 
decay again. And Nix's semimajor axis increases with 
time. 



2.3. Secular Evolution With Tidal Damping 
To explain the behavior of the simplified model. 



we 



consider first the effect of only the secular term [H = 
HsccCn), together with tidal damping (eq. [11]). The 
equations of motion are then 



d^ 
dt 



zc 
zn 



IT] 



S -(3S 
-0 1 



zc 
Zn 



ic 



Zc 




(13) 
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TABLE 1 

Terms in Charon-Nix Hamiltonian Used in This Paper'' 



2ac 

— —Mm. 

3unp,Af — 2a„ 



^^unp,C 

H2:1,CN = -IJ-c^ ig3Z*c + g4Z*^) e<(2^'V-Ac) +C.C. 

Hs-.i.CN = -MC^ {95Z*J + gez*^z*^ + g-jz*^) e^i^^^-Xc) + c.c. 



Charon's unperturbed energy 
Nix's unperturbed energy 
secular interaction energy 
2:1 resonant interaction energy 

3:1 resonant interaction energy 



"c.c." denotes complex conjugate of preceeding term. Coefficients gi are listed in Table 2. 



TABLE 2 

Coefficients in Ciiaron-Nix Hamiltonian'* 



-"''3/2/4 



91 

92 = 

53 = (-2 - aD/2)bj^.^ 
gi = (3/2 + ai)/2)fei/2 " 2a 

96 = (21/8 + 5aD/4 + a'^D'^ /9,)h\^^ 
gg = (-5 - 5aD/2 - a2£|2/4)52^^ 

97 = (17/8 + 5aD/4 + a"^ D"^ /9,)b\^.^ 
The coefficients arc given in Appendix B of Murray 

term for Jacobi coordinates (eq. [A23]). The Laplace coefficients 
numerical expressions, we set a = 4~^/^. 



- 27o/8 



: 0.0815 
: -0.0792 
: -0.390 

: 0.0811 

: 0.315 

: -1.41 

: 0.186 



Dcrmott (1999), except that we include the effective indirect 
are functions of a = ac/a-N, and D = d/da. In the 



where 

Mn y/a^ _„ 0-92 f... 

Mc \/ac 2.91 

Equation (13) is a linear equation with constant coeffi- 
cients. It is a simple exercise in linear algebra to solve 
it (Wu & Goldreich 2002). Since jc is much smaller 
than any of the frequencies in the problem, we first con- 
sider what happens in the absence of damping by setting 
7c = 0. Setting {zc,zn) oc e*"*, the eigenvalues and 
eigenf unctions are 




where 



(15) 
(16) 

(17) 



With non-zero 7c, both the + and — modes are damped. 
The damping rates are found by calculating the imag- 
inary part of the eigenvalues of the full ecjuation (eq. 
[13]), and expanding to first order in 7c, which yields 
the damping rates 



7± 



7c 
2 



(l±(l + e)" 



1/2 



(18) 



The physical meaning of the above results becomes 
clearer in the limit that 5 <C 1. For the eigenvalues 
and eigenf unctions in the absence of damping, consider 
first what happens when Charon's eccentricity is held 
fixed, zc =constant= zq- Then equation (13) shows 
that = fSzc + ZNjr, where l3zc is the constant forced 
eccentricity, and the free eccentricity zat,/,- has constant 
amplitude and a phase that precesses at frequency rj. 
Similarly, if one artificially held Nix's eccentricity fixed, 
Zn =constant= z^, then Charon would have forced ec- 
centricity = I3zn, and its free eccentricity would pre- 
cess at frequency rjS. When both zc and zn are allowed 



to evolve freely, there are two modes. In one of them, 
which can be called the CfN mode ("Charon forces Nix"), 
Nix's eccentricity is equal to its value forced by Charon, 
Zn = Pzc- From Charon's equation of motion, we see 
that this mode has frequency r]d{l — p-^), comparable to 
Charon's free precession frequency rid; therefore the CfN 
mode has 



zn 



oc 



/3> 



(19) 
(20) 



in agreement with the small S limit of the -|- mode in 
equations (15)-(16). 
In the other mode (NfC ="Nix forces Charon"), Nix 

has a free eccentricity, so it is freely precessing. Nix's ec- 
centricity tends to drive zc to precess around its forced 
value (3zn- But since Nix precesses much faster than 
Charon — by the factor 1/6 ^ 1 — Charon can only reach 
an eccentricity of \zc\ ~ (36zn in the time that Nix un- 
dergoes a full precession period. Since \zc\ <C \zn\, Nix's 
equation of motion is dzj^/dt « irjZN, showing that this 
mode has frequency equal to Nix's free precession fre- 
quency 77. And Charon's equation of motion shows that 
Zc — —(35zn, so the NfC mode has 

W_=T? (21) 

(22) 



Zn 



-135 
1 



in agreement again with equations (15)-(16) for (5 <C 1. 

When tidal damping is active, the damping rate of the 
CfN mode in the 5 <C 1 limit is just that of Charon in 
isolation 

7+ = 7C • (23) 

For the NfC mode, we take advantage of the fact that 

\zc \ ^ \zn\- Nix's approximate equation dzj^ /dt « irfz^ 
implies zn = ^atb*''*, where zn is (nearly) constant. 
Charon's equation is then dzc/dt w —ir]P6zNe^'^*—jczc\ 
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its solution after a few damping times of the CfN mode 
{t > l/7c) is zc = + i^c /rDzNd"^*' ■ We now 

consider the full equation for zjv, and substitute into 
this equation both the above zc and the zeroth order 
solution zn = ZnB^'^*, where zn is now allowed to be 
time-dependent, yielding dz^/dt — {irjfi'^S — jcSP'^)zn- 
Therefore \zn\ damps at the rate 

7- = IcSp^ . , (24) 
in agreement with the small 6 limit of equation (18). This 
gives the damping time quoted in the introduction (eq. 
[I])- 

To summarize, the purely secular evolution has two 
normal modes, one of which (CfN) decays very quickly. 
The second mode (NfC) is much more slowly damped. 
As shown in the introduction, the timescale to damp this 
mode can be shorter than the age of the Solar System for 
reasonable values of 7c and Mn/Mc- 

2.4. Adding in the 3:1 Resonant Forcing Term 

The purely secular evolution changes markedly when 
the 3:1 forcing term H^-\^cn is included. In Appendix 
B, we show that including the 3:1 term between Charon 
and Nix alters the damping rate of the slowly damped 
NfC mode from 7_ (eq. [24]) to 73:1 (eq. [B14]). Figure 
3 plots the ratio "f^-.i/l- as a function of the semimajor 
axis ratio. Remarkably, the 3:1 term has an effect at large 
distances from the location of the nominal 3:1 resonance, 
{aM/acf/^ = 3. 

Far beyond the nominal 4:1 resonance, {un / ac)^^"^ ^ 
4, the 3:1 resonant forcing term has little effect 
(73:1/7- ^ !)• But at Nix's current location (aj^/ac — 
4^/"^), the damping rate is reduced to only 7% of the 
purely secular rate. Furthermore, in the range 3.4 < 
(ajv/ttc)^^^ < 3.9 the damping rate is negative: instead 
of damping there is exponential growth. For still smaller 
values of (inIo-c the damping rate is very large and pos- 
itive. These behaviours are also confirmed by the full 
N-body integration (Fig. 1). 

2.5. Effect of the 2:1 Resonant Forcing Term 

In this subsection, we solve the simplified model of §2.2 
when Charon and Nix are coupled both secularly and via 
the 2:1 resonance, and tides damp Charon's eccentricity. 
We discard the 3:1 resonant forcing term. The equations 
for Ac and A at are given by equation (9), the eccentricity 
equations are given by equation (13) with the following 
term added to the right-hand side. 



d 
dt 



Zc 
Zn 



2:1 



= IjJLcnN 



53 <5 
94. 



J(2\n-Xc) 



(25) 



and the semimajor axis equations are 



d_ 
dt 



ac 



-Sac 
2ajv 



{93Z*c + 54^^)e'(2A«-Ao) + c.c. (26) 



Since ac,aN vary by only a small amount {O^ficCN) 
or smaller) on the 2:1 timescale, we may treat ncnpf 
as constants in the equations for the mean longitudes, 
which are then trivially solved 

Ac = net + const (27) 

Aat =njvi + const . (28) 




Fic;. 3. — Effect of the 3:1 resonance on the eccentricity damp- 
ing rate: The curve sfiows tfie ratio of 73;! (tlic secular damp- 
ing rate including the 3:1 resonance) to 7— (the rate without 
the 3:1 resonance), where 73;! is given by equation (B14) and 
7— by equation (24). In the range 2.3 < ajv/ac < 2.5 (i.e., 
3.4 < (ajv/ac)'*^^ < 3.9) the damping rate is negative, implying 
exponential growth. The x's show results from numerical integra- 
tions of equation (Bl), showing good agreement with the analytic 
expression. 

The general solution of the eccentricity equation is the 

sum of the two homogeneous (or "free") solutions given 
in §2.3 and the particular (or "forced") solution that is 
driven at frequency 



722:1 = 2nAr - nc ■ 



(29) 



Discarding the rapidly damped free solution (CfN), 
leaves 



Zc 
Zn 



= ZN,fr{t) 
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+ 



Zc 
Zn 



2:1 



where ZN,fr{t)=constxe'^'^-'i''^'^^^* and 



^C,2:l - 



--l^cg^o e 

^2:1 



in2:it 



l + i 



IC 
"2:1 



^iV,2:l = MC54 — e*"-i* 
"2:1 



(30) 

(31) 
(32) 



to first order in 7c. Even though the 0(7c) correction to 
zc,2:i is very small in absolute value, it plays an impor- 
tant role in the evolution of the semimajor axes because 
of its imaginary coefficient. Secular terms have been dis- 
carded from the forced solution, which is appropriate be- 
cause r] <C |n2:i|. 

We may now substitute these eccentricities into equa- 
tion (26). The free eccentricities produce rapidly oscil- 
lating terms oc e*("2^i+''^* that lead to small variations of 
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Charon's tidal damping time tc = l/7c- This is too 
short to be seen in Figure 1. The second normal mode 
(NfC=Nix forces Charon) is damped on a much longer 
timescale, at the rate 7_ = "fcSp"^ ^ 'yc^N /Mc- In 
this mode, Charon's proper eccentricity is forced by Nix 
to 6(7 ^ {Mj^/Mc)eN- But purely secular effects do not 
suffice to explain the evolution seen in Figure 1. The 2:1 
forcing term causes Nix's semimajor axis to increase at 
the rate k ~ 0.0067_. The 3:1 forcing term changes the 
damping rate of the NfC mode. At early times in Fig- 
ure 1, the 3:1 term forces the damping rate to be much 
faster than 7_. At later times, as Nix's semimajor axis 
increases, it enters a region where the 3:1 term causes the 
eccentricities to grow, rather than damp (see Fig. 3). At 
even later times, the proper eccentricities decay again, 
though at a rate less than 7_. 

Although we have ignored Hydra in our discussion, we 
describe its dynamics in Appendix C. 

3. DISCUSSION 
3.1. Current Orbits of the Moons 

Can Nix's current orbit be explained as the end state of 

the evolution seen in Figure 1? To answer this question, 
we first need a better understanding of the current orbits 
of the moons. 



Fig. 4. — Numerical simulation of Charon and Nix, in which 
Charon and Nix arc coupled both secularly and via the 2:1 res- 
onance, and tides damp Charon's eccentricity. Top panel shows 
eccentricities (ejv = \zn\,&g = \zc\) and bottom panel shows 
semimajor axes relative to Charon's initial semimajor axis ac o, 
confirming the solution given in equations (30) and (33). In the 
top panel, the slowly damped mode damps at the rate 'yc&0^-, 
and after 2-3 damping times the eccentricities reach their reso- 
nantly forced values; the dashed lines labelled |2;jv,2:i|, |2C,2:i| are 
the forced eccentricities given by equations (31)-(32). The dashed 
line in the bottom panel shows that Nix's semimajor axis increases 
at the rate k (eq. [34]). The addition of the 3:1 resonant forcing 
term will qualitatively change the eccentricity evolution, but not 
the semi-major axis evolution. 

ac,aN- But the forced eccentricities have a more inter- 
esting effect: they induce a slow secular change, 

dUc\J-Sac/2\ 
dt \aN J \ aN J 

K = 4:S^c ifJ-cgs—] ■ (34) 

In the time it takes Nix's free eccentricity to damp, 

its semimajor axis increases by the factor k/{0^5^c) = 
0.6/i^ ~ 0.6% at ajv = 4^/'^oc; in continues to increase 
even after the free eccentricities have damped away. 

Figure 4 shows a numerical integration of the equa- 
tions of motion given in this subsection, with the same 
parameters for the N-body integration shown Figure 1. 
The evolution seen in the Figure confirms the solutions 
given in equations (30) and (33). 

2.6. Summary of Pluto, Charon, and Nix's Evolution 

We have explained the behavior seen in the N-body 
simulation of Figure 1. There are three important types 
of interactions between Charon and Nix: secular, 3:1 
forcing, and 2:1 forcing. Secular interactions lead to 
two damped normal modes. One of these (CfN=Charon 
forces Nix) is rapidly damped away on the timescale of 



3.1.1. Nix's Eccentricity 

Buie et al. (2006) fit Nix's orbit to a Keplerian orbit 
and found ejv = 0.0023(21). Tholen et al. (2007b) fit all 
the moons' orbits simultaneously to the results of 4-body 
integrations, and found for the best-fit solution that Nix's 
eccentricity varied in time, with < cat < 0.0272. 

As seen in Figure 1, the instantaneous values of Nix's 
eccentricity (shown as points), can be much larger than 
the time-averaged eccentricity vector (shown as a line, 
for an averaging time of 100 days). A clearer view of 
Nix's orbital state at the end of that simulation may 
be seen by taking a Fourier transform of the eccentric- 
ity vectors around this time (bottom panel of Figure 5). 
The forest of high-frequency peaks at w > 0.1/day are 
forced eccentricities. They are associated with the res- 
onant terms in the disturbing function. For example, 
Nix's peak at w = \2nN — nc\ is the 2:1 forced eccentric- 
ity described above (eq. [32]). The low- frequency peak at 
ui ~ 0.0035/day is Nix's secular eccentricity (which may 
also be called its "free" or "proper" eccentricity). Even 
though the forced eccentricities dominate the instanta- 
neous eccentricity, they are not relevant if one wishes to 
use Nix's current orbital state to infer something about 
its past. This is because for fixed semimajor axes and 
masses, the forced eccentricities are fixed. It is only the 
secular eccentricity vector that represents a true degree 
of freedom. Tides act to damp the secular eccentricity, 
but they have no effect on the forced eccentricity. At 
the time depicted in Figure 5 the CfN mode has damped 
away, and only the NfC mode remains. The frequency of 
the NfC mode is U- = rj (eq. [21]), or w_ - 0.004/day 
in agreement with the secular peak mxii in the figure, i.e., 
the low-frequency peak is the remains of the NfC mode 
as it is being damped by tides. 

In Figure 6, we "observe" the orbital state at the end 
of the simulation in the manner of Buie et al. (2006) by 
fitting Nix's orbit to a Keplerian ellipse around a point 
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Fig. 5. — Fourier decomposition of Charon & Nix's eccentric- 
ity vector for the final state in Fig. 1. The Fourier amplitude 
= J {ecoszu) cos{Lijt)dt with orbital elements measured in Ja- 
cobi coordinates. The eccentricity of Charon is dominated by the 
2:1 forcing term, while Nix by a combination of short term forcing. 
The peaks (indicated by aj„) are associated with the slowly decay- 
ing secular mode (NfC) which precesses at aj_ (eq. [22]). Their 
amplitude is a factor of ~ 100 below the total eccentricity. 

mass at the barycenter with mass Mp + Mc (where the 
value of Mp + Mq is to be found by the fit). We also 
fit Charon's orbit relative to Pluto with a Keplerian el- 
lipse. We output the positions of Charon and Nix as 
functions of time from our numerical integrator. We then 
use the downhill simplex method (Press et al. 1992) to 
search fits for the set of 5 parameters: ecostu, esinw, 
a, Mp + Mc and r (epoch of periapse passage). We find 
that the Keplerian-fit eccentricities are the secular parts 
of the total eccentricity, with the instantaneous values 
being much greater but averaged out in the fit. There- 
fore, somewhat counterintuitively, the Keplerian fit of 
Buie et al. (2006) provides a more useful diagnostic of 
Nix's eccentricity than the 4-body fit of Tholen et al. 
(2007b). 

3.1.2. Charon's eccentricity 

Tholen et al. (2007b) report that ec = 0.00348(4) 
when they do either a 4-body or a Keplerian fit, and 
that their earlier, much smaller value for ec (Buie et al. 
2006) was incorrect. We find this large value of ec very 
puzzling. Our theory predicts that Charon's eccentric- 
ity rapidly damps away by tides. More precisely, on the 
timescale rc < 10 Myr, the CfN mode damps away. Af- 
ter this happens, Charon's secular eccentricity in the NfC 
mode is ec ~ {Mm /Mc)eM as seen also in Figures 5-6. 
Thus Charon's eccentricity should be much smaller than 



Fig. 6. — Results of Charon & Nix orbital fitting for the final state 
in Fig. 1, insisting on Keplerian ellipses . We output data (points 
for Charon and open circles for Nix) once a month for 10 years, 
and do a Keplerian fit every 12 months (solid triangles for Charon 
and open triangles for Nix, error-bars in time-axis indicating span 
of Keplerian fit). Eccentricities are shown in the left panel: excp 
fall well below the instantaneous values and are approximately the 
secular values (obtained from Fig. 5). The upper-right panel shows 
the value of ro, which varies on orbital timescales and are not 
captured by the Keplerian fits. The lower-right panel displays the 
fitting residuals measured by r — r^op in unit of Rp. Charon's 
residuals are magnified by a factor of 100. 

the Tholen et al. (2007b) value. ^ It is highly implausi- 
ble that Charon's eccentricity was excited in the last 10 
Myr by some external event, such as fiybys by passing 
Kuiper belt objects (Stern et al. 2003). In our view, the 
most plausible resolution of this puzzle is that the high 
eccentricity of Tholen et al. (2007b) is incorrect; perhaps 
inhomogeneities on the surface of Pluto or Charon are 
responsible for giving an apparent eccentricity. Future 
observations should help to resolve the puzzle. 

3.2. Constraints on Nix's Initial Orbit 

Let us consider the following scenario for explaining 
why Nix currently has a very low proper eccentricity, 
eN ^ 0.002 (Buie et al. 2006): perhaps Nix formed with 
a high eccentricity at ^ 2.25ac, and then tidal evolu- 
tion circularized its orbit as it was migrated to its current 
location, as in Figure 1. Note that Nix's proper eccentric- 
ity at the end of that simulation is close to the Keplerian- 
fit value of Buie et al. (2006) . The bottom panel of Figure 

^ Our theory predicts that Charon's secular eccentricity should 
be much too small to be observable. But its forced eccentricity 
should be ~ 0.3Afjv/A/p ~ 10~^, forced primarily by the 2:1 reso- 
nance with Nix, and hence rapidly precessing at the 2:1 frequency 
(eq. [31], Figs. 5-6). If this forced eccentricity is measured, one 
could infer from it Nix's mass. 
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Fig. 7. — Trajectory' of Nix's orbit when tides damp Charon's 
eccentricity. The top panel shows a numerical integration of the 
simplified model. The simulation is identical to the one shown in 
Figure 2, but here the eccentricity is plotted versus semimajor axis 
ratio. Similarly, the solid line in the bottom panel is the proper ec- 
centricity from the N-body simulation shown in Figure 1, with the 
value of ajv averaged over 100 days. The other curves in the lower 
panel are from N-body simulations that started with Nix both fur- 
ther in and further out, but also with ejv = 0.05. The simplified 
model captures the essential behavior of the system but differs in 
detail, partly because the Laplace coefficients arc taken to be con- 
stants, and partly because our disturbing function is only accurate 
to 0(iJ,c)- Nix's current proper eccentricity is 0.002 (Buic ct al. 
2006). The nominal location of the 4:1 resonance is shown as an 
arrow. Trajectories that have very low eN,soc do not experience 
an appreciable jump when crossing this resonance (confirmed by 
further integration not shown here), while trajectories that have 
higher Csec (> 0.02) do. 



7 shows two additional N-body simulations with differ- 
ing initial ajv- It illustrates that if Nix started inward 
of ~ 2.25ac, its proper eccentricity is damped when it 
reaches its current location near the 4:1 resonance. But 
if it started beyond 2.2hac-i it would have retained much 
of its initial eccentricity. 

The difficulty with this scenario is that to have mi- 
grated Nix from inward of 2.25ac to its current position, 
we require that tc < 2 x 10Vrs(/ijv/0.00015) (Fig. 1), 
where /xjv = M^/Mp. Even for /ijv ^ 0.00015, this re- 
quires a tidal damping > 10 times more efficient than the 
current estimate for the Phito-Charon binary (Dobrovol- 
skis et al. 1997). This seems difficult unless Charon is 
a molten sphere with large Love number k2 (tidal dis- 
tortion) and small Q value. Moreover, it is a strange 
coincidence that Nix is so close to the 4:1 location today, 
since in our theory the outward migration of Nix does 
not pause at 4:1. 

So unless the tidal dissipation time of Charon tc % 10^ 
yrs. Nix will have to be initially deposited at its present 
semimajor axis with proper eccentricity 0.002, its current 
measured value. This places very stringent constraint for 
the formation scenario and rules out all but formation in 
a disc. In the latter case, the near 4: 1 location is a result 
of migration in the disc (Shannon et al., in preparation). 

Conversely, the fact that Nix could not have migrated 
further than from inward of the inner-most stable orbits 
(ajv ~ 2.15ac) restricts rc to be > 2 x 10^ yrs. 



APPENDIX 

APPENDIX A. PLANETARY EQUATIONS IN JACOBI COORDINATES 
The Hamiltonian for Pluto, Charon, and Nix is 

P| P% _ GMpMc _ GMnMc _ GMmMp 

2Mp ^ 2Mc ^ 2Mn ~ \Rp - Rc\ ~ \Rn-Rc\~ \Rn - Rp\ 



H ■ 



(Al) 



where Rj are position vectors from an arbitrary inertial origin, Pj = MjdRj/dt are the conjugate momenta, and Mj 
the masses. In Jacobi coordinates, Charon's position is measured relative to Pluto's, and Nix's position is measured 
relative to the center-of-mass of the Charon-Pluto binary. We transform to Jacobi coordinates in two steps. First, we 
employ the generating function 



F={Rc- Rp)-Ppc + 



McRc + MpRp 
Mc + Mp 



PC 



(A2) 



to switch the coordinates of the Pluto-Charon binary from Rp, Rc, Pp, Pc to the binary's center-of-mass and its 
relative position vector. 



R 



PC- 



McRc + MpRp 



Mc + Mp 
rpc = Rc-Rp , 



(A3) 
(A4) 
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and their conjugate momenta {Ppc,Ppc)i yielding the new Hamiltonian 

Ppc , Ppc , Pn G{Mp + Mc)Mpc GM^Mc GM^Mp 

= , ^ — : — , , X + TTTT h 



2{Mp + Mc) 2Mpc 2Mn rpc \Rn - Rpc - - IJ-c^pcI \Rn - Rpc + IJ-cr pc 
where 

= (A6) 

is the reduced mass and 

Mc 

= Mp + Mc • ^^^^ 

Prom Hamilton's equation for dRpc/dt, drpc/dt, we see that Ppc is the total momentum of the binary, and ppQ is 
the momentum of the relative orbit. 
We complete the transformation to Jacobi coordinates with the generating function 

F = Rpc-PpcN + {Rn — Rpc)'Pn i-^^) 

which yields the new coordinate 

rN = RN- Rpc , (A9) 

as desired, as well as the new coordinate RpcN = Rpc\ the conjugate momenta are pj^ = Pn and the total 
center-of-mass momentum 

PpcN = Ppc + Pn ■ (AlO) 

Since this is constant {Rpcn does not appear in the Hamiltonian), we may set it to zero. The Hamiltonian in Jacobi 
coordinates is 

H(.,...Pc;p„..„)=J5^ + 4(l + -^) (AH) 

G{Mp + Mc)Mpc GMnMc GMnMp 

rpc \rN - - fJ'c)rpc\ {tn + M-crpcl ' 

Thus far. our treatment has been exact. Henceforth, wc drop the factor Mm /{Mp + Mc) in the second term above; 
although it is simple to retain it, we drop it for notational convenience. We choose units so that 

G{Mp + Mc) = 1 . (A13) 

By adding and subracting the term —Mn/tn, the Hamiltonian may be written as 

H{Ppc, rpc; Pn^ '^n) = H^np,PC + H^np,N + Hpert , (A14) 



where 



rr _ Ppc Mpc , . 

•Hunp.PC - 7^ (A15) 

ZMpc rpc 

_ p% Mn , . 

MnMc ( 1 l\ MnMp ( 1 1\ 



Mp + Mc \\rN - - fJ-c)rpc\ tn J Mp + Mc \\rN + ^^crpc\ tn 

Since we seek equations for the orbital elements, we transform variables to {apc-i)^pci^pc-,''^pc} and 
{aAT, Aat, Cat, Wat}, defined in the usual way. (The mass of the cental body that enters into the definitions is 
Mp + Mc, both for Nix and for Pluto- Charon.) The equations of motion Lagrange's planetary equations 
are given by Hamilton's equations for canonical variables, which we may take to be the Poincare variables , e.g. 
Kpc = Mpc \/a.pc , A AT = Mn^/on, etc. (Murray & Dermott 1999). It remains to express the Hamiltonian in terms 
of the orbital elements. The unperturbed terms are 

/?unp,PC = -^ (A18) 

Hunp.N = -^ (A19) 

For the perturbed term, we resort to the usual Fourier expansion into a sum of cosine terms. Appendix B of Murray 
& Dermott (1999) tabulates the coefficients for the Fourier expansion of l/|r — r'\ in terms of the orbital elements of 
one body at position r and an exterior body at position r'. (More precisely, that Appendix tabulates the direct part 
of the disturbing function, TZd = a' /\r — r'\.) Since equation (A17) consists of four terms of this form, we can easily 
extract the coefficients from Murray & Dermott (1999). These coefficients are functions of a, the ratio of semi-major 
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axes of the body at r to the one at r'. The coefficient for the term l/\rN — (1 — Mc)'"pc| should be evaluated at 
a = (1 - nc)apc,N, where 

apc,N = — ; (A20) 

similarly, the coefficient I/tn should be evaluated at a = 0, and that for l/\rN + ficfpcl at a = hcC(pc,n- But there 
is an extra complication with the latter term, since fjbcfpc enters with a positive sign, instead of a negative one. To 
correct for this, one must set wpc — > t^pc + and Apc ~* •^pc + in the argument of the cosine; equivalently, one 
must multiply by —1 if within the cosine argument the sum of the integer coefficients of vopc and of Xpc give an odd 
number. 

The procedure described above yields the equations of motion as precisely as desired when enough Fourier terms 
are retained in the disturbing function iJpert (aside from the factor Mjv/(Mp + Mc) dropped from the Hamiltonian.) 
But for the purposes of this paper, it suffices to obtain coefficients that are incorrect by ^ /ic ^ 10%. Therefore we 
expand Hp^rt to leading order in fic- For each cosine term of the direct potential TZd = aN/\rN — rpc\ of the form 

1lD{apc,N) cos (j) (A21) 

(suppressing the other arguments ofTZD), we have 

^^pert = -Hc—Tlcos (f) , (A22) 

where 

n = noiapcN) - TZoiO) ± apc,NiDTZD)\a=o , (A23) 
with D = d/da and ± — > + if the integer d(p/dXpc + d<p/dzupc is even; otherwise ± — > — . 

APPENDIX B. EFFECT OF THE 3:1 RESONANCE ON THE SECULAR DAMPING RATE 

We solve the equations of motion for zc and zn, including secular interactions, the 3:1 resonance, as well as tidal 
damping. The equations are given in §2.2, except here we discard the term H2:i,cn', explicitly, 

i(-)--<(_vfO(s)-K-r?oci)'""--^"--Co") • 

where 

v = 2nNlicg7, X=-p-, v= — , (B2) 

and the other quantities are the same as in equation (13). To solve these equations analj^tically, we take the semimajor 
axes to be fixed, setting 

3Ajv - Ac = n3;it (B3) 

where 

713:1 = 3njv - nc (B4) 

is constant. We solve the equations perturbatively in the small parameter (5 <C 1, adopting the same approach as in 
the paragraph above equation (24), i.e., we assume that \zc\ "C |^jv|, which may be verified a postiori.® To leading 
order in \zc\, 

^ ~ irjZN + ivz*t,e'^'-'* . (B5) 
Replacing the ~ with =, this equation has solution 

= fce*(''+';)* + ^fc*e*("-^^i-''-'')* , (B6) 

V 

where k is the complex integration constant, and b is either of the two roots of the quadratic equation that results 
from 

^ = — V—r. ■ m 

V n3:i — 277 — 6 
For definiteness, we choose the low-frequency root 




(B8) 

In the vicinity of Nix's current semimajor axis n^-.i < 0, implying that b ~ /n^-.i to leading order in \ri/nz:i \ ~ 2% 
and |z^/n3:i| 4%. 

* If initially \zc\ ~ |^iv|, then zc would quickly decay on timescale to a value ~ i5|zjv|. 
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Inserting equation (B6) into the approximate equation for zc, 

dzc 



dt 



yields 



Zc = -5-. 



+ T]-ijc 



"3:1 -b-T]-ijc 
-1 



_l^*g,i{n3:i-b-ri)t 



(B9) 



(BIO) 



discarding the homogeneous solution since it decays away after time ~ 

Next, we rewrite the full equation for zn by inserting the approximate solution (B6) into equation (Bl), with k now 
time-varying, resulting in 

rik p-i(b+v)t 

di = 62/^2 [ivP + bx)zc + ivWiy + iyx)zhe'''^-^'] (Bll) 

Upon substitution of equation (BIO), we arrive at 

dk 

where 

5 

V 



ipk + iqk^e'^''^-'-^^-^"^* 

{nP + hxf 



+ 



{■n(3b/v + vxY 



b + rj-i^c n3;i - b - T] + ijc 



(B12) 



(B13) 



and (J is a complex constant whose explicit form we do not give because we shall have no use for it. Equation (B12) 
has the same form as equation (B5), except that now the coefficients are complex. Its solution has the same form as 
equation (B6), except now multiplied by the prefactor e~P™(^^l*. In conclusion, on timescales much longer than 
Zn slowly decays at the rate 73:1 = Im(p), or 



73:1 



'c 



V/3 + bx 
b + rj 



■I] 0b/ V + vx 

ri3:l -b-T] 



(B14) 



1-62/1/2 

This rate is plotted in Figure 3 relative to damping rate in the absence of the 3:1 forcing term. 

APPENDIX C: HYDRA 

Hydra's evolution in the presence of Chaon is similar to Nix's, but with some quantitative differences. Hydra's 
purely secular damping rate relative to Nix's is (eq. [18]) 



1-N 



M 



N 



0.7 



N 



■Ml 

M. 



(CI) 



where subscripts N and H arc for Nix and Hydra, and the numerical expression is evaluated for Nix and Hydra at, 
respectively, the nominal 4:1 and 6:1 resonances with Charon. The 3:1 resonance has a relatively small effect on 
Hydra's secular damping rate: from the extrapolation of Figure 3 to the nominal 6:1 resonance, it reduces 7_ij by the 
factor 0.6. Hydra's forced eccentricity due to the 2:1 resonance with Charon is smaller than Nix's by (eq. [32]) 



ZH,2-1 



94:,H riH n2:l,N 
94,N nN n2:l,H 



0.2 



and its migration rate is reduced by (eq. [34]) 

Kh_ _ MHy/aH ( 93,H riH n2:l,N 

kn Ma 



Ny'dN \93,N n2:l,H 



0.09 



■Ml 

M, 



H 



N 



(C2) 



(C3) 



It is also interesting to consider the interaction between Nix and Hydra. Hydra's reported eccentricity based on a 
Kcplcrian fit differs from zero, en = 0.0052(11), and this might be due to the proximity of Nix and Hydra to their 
mutual 3:2 resonance. The 3:2 interaction energy between Nix and Hydra is 



H3..2,NH = -MiV^ {9aZ% + 9tZ*„) e'(3AK-2A„) ^ ^ ^ 



(C4) 



2/3 



where hn = Mn/{Mc + Mp), 9a = -(6 + aD)bl^^/2 = -2.03, and 5^ = (5 + aD)bl^^/2 = 2.48, using a = (2/3) 
in the numerical expressions. Using the equation for Hydra's eccentricity, it is simple to show that the 3:2 resonant 
term gives a forced eccentricity to Hydra equal to 



Zh,3:2 = IJ-Ngt 



riH 



J(3\h-2Xn) 



(C5) 
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With the observationally derived values for jih and tin given in the Introduction, we find |zif,3:2| = 0.003/iAr/lO^^, 
close to Hydra's reported eccentricity. Since Nix and Hydra are close to 3:2 resonance, this eccentricity vector precesses 
at a slower rate than other resonantly forced eccentricity vectors (which precess at orbital timescales). It will show up 
in a Keplerian fit that covers many orbital periods. And we suggest that this likely explains the observed eccentricity 
of Hydra. Nix's 3:2 forced eccentricity has a similar expression but scales with fin. Tholen et al. (2007b) report a 
Hydra mass that is lower by a factor of 2 than Nix and encompasses zero. This is consistent with Nix's lower measured 
eccentricity. 
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